home
***
CD-ROM
|
disk
|
FTP
|
other
***
search
/
Fritz: All Fritz
/
All Fritz.zip
/
All Fritz
/
FILES
/
PROGMISC
/
PCSSP.LZH
/
PC-SSP.ZIP
/
MATINSL.ZIP
/
GELS.FOR
< prev
next >
Wrap
Text File
|
1985-11-29
|
5KB
|
183 lines
C
C ..................................................................
C
C SUBROUTINE GELS
C
C PURPOSE
C TO SOLVE A SYSTEM OF SIMULTANEOUS LINEAR EQUATIONS WITH
C SYMMETRIC COEFFICIENT MATRIX UPPER TRIANGULAR PART OF WHICH
C IS ASSUMED TO BE STORED COLUMNWISE.
C
C USAGE
C CALL GELS(R,A,M,N,EPS,IER,AUX)
C
C DESCRIPTION OF PARAMETERS
C R - M BY N RIGHT HAND SIDE MATRIX. (DESTROYED)
C ON RETURN R CONTAINS THE SOLUTION OF THE EQUATIONS.
C A - UPPER TRIANGULAR PART OF THE SYMMETRIC
C M BY M COEFFICIENT MATRIX. (DESTROYED)
C M - THE NUMBER OF EQUATIONS IN THE SYSTEM.
C N - THE NUMBER OF RIGHT HAND SIDE VECTORS.
C EPS - AN INPUT CONSTANT WHICH IS USED AS RELATIVE
C TOLERANCE FOR TEST ON LOSS OF SIGNIFICANCE.
C IER - RESULTING ERROR PARAMETER CODED AS FOLLOWS
C IER=0 - NO ERROR,
C IER=-1 - NO RESULT BECAUSE OF M LESS THAN 1 OR
C PIVOT ELEMENT AT ANY ELIMINATION STEP
C EQUAL TO 0,
C IER=K - WARNING DUE TO POSSIBLE LOSS OF SIGNIFI-
C CANCE INDICATED AT ELIMINATION STEP K+1,
C WHERE PIVOT ELEMENT WAS LESS THAN OR
C EQUAL TO THE INTERNAL TOLERANCE EPS TIMES
C ABSOLUTELY GREATEST MAIN DIAGONAL
C ELEMENT OF MATRIX A.
C AUX - AN AUXILIARY STORAGE ARRAY WITH DIMENSION M-1.
C
C REMARKS
C UPPER TRIANGULAR PART OF MATRIX A IS ASSUMED TO BE STORED
C COLUMNWISE IN M*(M+1)/2 SUCCESSIVE STORAGE LOCATIONS, RIGHT
C HAND SIDE MATRIX R COLUMNWISE IN N*M SUCCESSIVE STORAGE
C LOCATIONS. ON RETURN SOLUTION MATRIX R IS STORED COLUMNWISE
C TOO.
C THE PROCEDURE GIVES RESULTS IF THE NUMBER OF EQUATIONS M IS
C GREATER THAN 0 AND PIVOT ELEMENTS AT ALL ELIMINATION STEPS
C ARE DIFFERENT FROM 0. HOWEVER WARNING IER=K - IF GIVEN -
C INDICATES POSSIBLE LOSS OF SIGNIFICANCE. IN CASE OF A WELL
C SCALED MATRIX A AND APPROPRIATE TOLERANCE EPS, IER=K MAY BE
C INTERPRETED THAT MATRIX A HAS THE RANK K. NO WARNING IS
C GIVEN IN CASE M=1.
C ERROR PARAMETER IER=-1 DOES NOT NECESSARILY MEAN THAT
C MATRIX A IS SINGULAR, AS ONLY MAIN DIAGONAL ELEMENTS
C ARE USED AS PIVOT ELEMENTS. POSSIBLY SUBROUTINE GELG (WHICH
C WORKS WITH TOTAL PIVOTING) WOULD BE ABLE TO FIND A SOLUTION.
C
C SUBROUTINES AND FUNCTION SUBPROGRAMS REQUIRED
C NONE
C
C METHOD
C SOLUTION IS DONE BY MEANS OF GAUSS-ELIMINATION WITH
C PIVOTING IN MAIN DIAGONAL, IN ORDER TO PRESERVE
C SYMMETRY IN REMAINING COEFFICIENT MATRICES.
C
C ..................................................................
C
SUBROUTINE GELS(R,A,M,N,EPS,IER,AUX)
C
C
DIMENSION A(1),R(1),AUX(1)
IF(M)24,24,1
C
C SEARCH FOR GREATEST MAIN DIAGONAL ELEMENT
1 IER=0
PIV=0.
L=0
DO 3 K=1,M
L=L+K
TB=ABS(A(L))
IF(TB-PIV)3,3,2
2 PIV=TB
I=L
J=K
3 CONTINUE
TOL=EPS*PIV
C MAIN DIAGONAL ELEMENT A(I)=A(J,J) IS FIRST PIVOT ELEMENT.
C PIV CONTAINS THE ABSOLUTE VALUE OF A(I).
C
C
C START ELIMINATION LOOP
LST=0
NM=N*M
LEND=M-1
DO 18 K=1,M
C
C TEST ON USEFULNESS OF SYMMETRIC ALGORITHM
IF(PIV)24,24,4
4 IF(IER)7,5,7
5 IF(PIV-TOL)6,6,7
6 IER=K-1
7 LT=J-K
LST=LST+K
C
C PIVOT ROW REDUCTION AND ROW INTERCHANGE IN RIGHT HAND SIDE R
PIVI=1./A(I)
DO 8 L=K,NM,M
LL=L+LT
TB=PIVI*R(LL)
R(LL)=R(L)
8 R(L)=TB
C
C IS ELIMINATION TERMINATED
IF(K-M)9,19,19
C
C ROW AND COLUMN INTERCHANGE AND PIVOT ROW REDUCTION IN MATRIX A.
C ELEMENTS OF PIVOT COLUMN ARE SAVED IN AUXILIARY VECTOR AUX.
9 LR=LST+(LT*(K+J-1))/2
LL=LR
L=LST
DO 14 II=K,LEND
L=L+II
LL=LL+1
IF(L-LR)12,10,11
10 A(LL)=A(LST)
TB=A(L)
GO TO 13
11 LL=L+LT
12 TB=A(LL)
A(LL)=A(L)
13 AUX(II)=TB
14 A(L)=PIVI*TB
C
C SAVE COLUMN INTERCHANGE INFORMATION
A(LST)=LT
C
C ELEMENT REDUCTION AND SEARCH FOR NEXT PIVOT
PIV=0.
LLST=LST
LT=0
DO 18 II=K,LEND
PIVI=-AUX(II)
LL=LLST
LT=LT+1
DO 15 LLD=II,LEND
LL=LL+LLD
L=LL+LT
15 A(L)=A(L)+PIVI*A(LL)
LLST=LLST+II
LR=LLST+LT
TB=ABS(A(LR))
IF(TB-PIV)17,17,16
16 PIV=TB
I=LR
J=II+1
17 DO 18 LR=K,NM,M
LL=LR+LT
18 R(LL)=R(LL)+PIVI*R(LR)
C END OF ELIMINATION LOOP
C
C
C BACK SUBSTITUTION AND BACK INTERCHANGE
19 IF(LEND)24,23,20
20 II=M
DO 22 I=2,M
LST=LST-II
II=II-1
L=A(LST)+.5
DO 22 J=II,NM,M
TB=R(J)
LL=J
K=LST
DO 21 LT=II,LEND
LL=LL+1
K=K+LT
21 TB=TB-A(K)*R(LL)
K=J+L
R(J)=R(K)
22 R(K)=TB
23 RETURN
C
C
C ERROR RETURN
24 IER=-1
RETURN
END